clear
clf

f = 400; b = 0.1;
P = 9; R = 0.09;
x_ = 20;

nSample = 1e6;
xTrue = normrnd(x_, sqrt(P), nSample, 1);
yMeas = normrnd(f*b ./ xTrue, sqrt(R), nSample, 1);

x = 10:.03:26;
J = 1/(2*R).*(yMeas - f * b ./ x).^2 + 1/(2*P).*(x_ - x).^2;
[JMin, xMAPId] = min(J, [], 2);
xMAP = x(xMAPId)';

histogram(xMAP, 'Normalization', 'pdf', 'BinWidth', .21)

hold on
plot([x_,x_],[0,0.25])
eXN = mean(xMAP);
plot([eXN, eXN], [0, 0.25])

eMean = eXN - x_
eMean = mean(xMAP - x_)

eSq = mean((xMAP - x_).^2)


